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Abstract. In this article, we study the application of Multi-Level Monte 
Carlo (MLMC) approaches to numerical random homogenization. Our ob- 
jective is to compute the expectation of some functionals of the homoge- 
nized coefficients, or of the homogenized solutions. This is accomplished 
within MLMC by considering different levels of representative volumes 
(RVE), and, when it comes to homogenized solutions, different levels of 
coarse-grid meshes. Many inexpensive computations with the smallest RVE 
size and the largest coarse mesh are combined with fewer expensive com- 
putations performed on larger RVEs and smaller coarse meshes. We show 
that, by carefully selecting the number of realizations at each level, we can 
achieve a speed-up in the computations in comparison to a standard Monte 
Carlo method. Numerical results are presented both for one-dimensional 
and two-dimensional test-cases. 

1. Introduction 

Many multi-scale problems have uncertainties at the smallest scales, that 
are due to the incomplete knowledge one has of the microstructure. For 
example, when considering porous materials, the microstructure is often 
generated based on some limited statistical information. This can lead 
to large uncertainties in terms of microscale heterogeneities. These un- 
certainties at micro-scales need to be mapped onto the simulations on a 
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coarse-grid, and this typically leads to considering large representative vol- 
umes (RVE) for these microstructures. 

In practice, the upscaled quantities that are used at the macroscopic 
level are computed using the solution of some local problems posed on 
these microstructures. It is often needed to solve many such local prob- 
lems (corresponding to many different random realizations, or snapshots, 
of the microstructure) , each of which being expensive due to the presence 
of small scales. The resulting amount of computational work may thus be 
prohibitively expensive. In this article, our objective is to design a com- 
putational approach that allows for fast calculations of the coarse-scale 
quantities based on fewer realizations. 

Our idea is to apply the Multi-Level Monte Carlo (MLMC) framework 
to multi-scale simulations. The MLMC approach was first introduced by 
Hcinrich in |23| for finite- and infinite-dimensional integration. Later on, it 
was applied to stochastic ODEs by Giles (see [HI [50]). More recently, this 
approach has been used for PDEs with stochastic coefficients by several 
authors, see [H [THl \H HH [57] . To compute an approximation of the expec- 
tation E(X) of some random variable X, the MLMC approach consists in 
considering several random variables Xi , at different levels I , that approx- 
imate X with various accuracies. The main idea is then to use different 
numbers of samples (i.e. independent realizations) at different levels. More 
precisely, many samples are used at the coarsest, less accurate level where 
the computation for each realization is inexpensive, while fewer samples are 
used at the finest, most accurate level that is expensive to compute. Com- 
bining the results of these computations by carefully selecting the number 
of realizations at each level can speed-up the computations in comparison 
to a standard Monte Carlo (MC) approach, where only one level (that of 
the quantity of interest itself) is considered. See Section l2~2l below for more 
details on the MLMC approach. 

In the framework of numerical stochastic homogenization, local problems 
are solved on representative volumes (RVE), and apparent effective prop- 
erties are next defined as averages of the solutions of these local problems 
over the RVEs. The computations on the RVEs are usually expensive, 
because large RVEs need to be considered to obtain effective properties 
with a reasonable accuracy. In the framework of MLMC approaches, our 
idea is to use RVEs of different sizes, and to consider many independent 
realizations of the smaller ones, for which the associated local problem is 
inexpensive to solve, and fewer realizations of the larger ones. 

The convergence of the MLMC approach depends on the accuracy of 
the computations at each level. Assessing how this accuracy improves 
when more expensive computations are considered is critical to determine 
how to choose the number of realizations at each level. In our case, we 
thus have to determine how the accuracy of apparent effective properties 
depend on the RVE size. Such estimations are not easy to obtain, both 
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from a theoretical and a practical viewpoint. In this work, we use the fact 
that, under some assumptions on the heterogeneous coefficients, it is known 
that the accuracy of the effective property approximation scales as (e/rf)P 
for some (3 > 0, where r\ is the RVE size and e is the characteristic small 
lengthscale of the heterogeneities (see e.g. [31 [5j |51 |TT 1 [T7 1 \W \ 13211231 135 ] ) . 

When the MLMC approach is used to compute the expectation of some 
functionals of the homogenized solution (rather than the homogenized co- 
efficient), we can use RVEs of different size to compute the homogenized 
coefficients, and also coarse grids with various size to solve the coarse scale 
equation. In addition to assessing the accuracy of the approximation of 
the effective properties in each RVE, we need to assess the accuracy when 
solving the coarse-scale equation. Standard FEM results are then useful. 

An important remark is that MLMC approaches are interesting when ef- 
fective properties are stochastic (otherwise, such approaches are as efficient 
as a standard MC approach). This situation appears in many applications, 
although homogenization theories for this case are less studied. Most ho- 
mogenization theories are indeed developed for ergodic coefficients that 
vary over a single scale. In this case, the apparent homogenized quantities, 
when computed on infinitely large RVEs, are deterministic. In the sequel, 
we briefly discuss homogenization results when the homogenized coefficient 
is stochastic (even when infinitely large RVEs are considered), and we use 
these results in our MLMC approach to adequately select the number of 
realizations at each level (namely, for each RVE size and each coarse grid 
size). 

Consider now the specific question of computing homogenized solutions 
with several grids of different size. For each of these grids, we first need 
to precompute the effective properties, say at each Gauss point of the 
macroscopic grid. Assume that these coarse grids are nested. Then, once 
the effective properties have been computed at the finest level (i.e. for the 
Gauss points of the finest grid) , no additional precomputation is needed to 
compute effective properties for the coarser grids (since their Gauss points 
are a subset of the Gauss points of the finest grid). In this case, we propose 
to use a weighted MLMC approach, where we give different weights to each 
level, so as to optimize the accuracy at a given cost. 

Our article is organized as follows. In Section [5J we briefly review theo- 
retical homogenization results and describe in details the MLMC approach 
in a general context. In Section [3J we next describe how to apply the 
MLMC approach to compute an approximation of the homogenized coef- 
ficients, and assess the accuracy of the proposed approach. We next turn 
in Section 3] to the computation of the homogenized solutions, using ei- 
ther the MLMC or the weighted MLMC approaches. Numerical results 
are collected in Section [5j We consider the case of the effective coefficients 
in Section 15.21 and of the homogenized solutions in Section 15.31 In both 
cases, we show that the MLMC approach yields a significant speed-up in 
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comparison to a standard MC approach. 



2. Preliminaries 

2.1. Numerical homogenization 

In this section, we describe the numerical homogenization procedure we 
use. Consider the problem 

- div(A e (x,uj)Wu e ) = f in D, (1) 

where D is an open bounded subset of M. d , A e (x,uj) is a heterogeneous 
random field (with a small characteristic length scale e), ui designates a 
random realization and / G L 2 (D) is a non-random function. We com- 
plement the problem (fTJ) with some boundary conditions that we do not 
specify, such that its solution u e is well defined (for instance, u e = on dD 
almost surely) . Furthermore, we assume that A e is uniformly bounded and 
coercive, in the sense that there exists two positive deterministic numbers 
< a m j n < a max such that, for any e, any £ € R d and any 1 < i,j < d, 
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almost everywhere in D and almost surely. 

For almost all realizations w, we consider a numerical homogenization 
procedure as follows. Given a representative volume centered at a macro- 
scopic point x with size 77, 

we solve, for any 1 < i < d, the local problems 

div(^(y,«)VXi(tf,w)) = Qiiil~ X2 (2/,^)=2/ l on5r t f. (2) 

Note that the precise boundary conditions used in these local problems 
are not essential when there is a scale separation. Rather than Dirichlet 
boundary conditions as in ([2]) , it is also possible to use Neumann boundary 
conditions, or periodic boundary conditions (see (TTT 125] ). 

Then, we define the apparent homogenized matrix A*(x,u)) by 

VI < i < d, A*(x,uj)e l = —7 [ A e (y,uj)Vxt(y,u) dy, 

T Jy* 

where is the unit vector in the direction % (i = 1, . . . , d). We denote this 
local homogenization procedure by 'H ri , i.e. 

A*(x,u) = H n {A e (x,uj)). 
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-div(^Vxi) = 




RVE 



Figure 1: Illustration of the numerical homogenization procedure. 



This procedure is repeated at every macroscopic point (see Figure [T] for 
illustration). Then, the coarse-scale equation associated to (UJ) is 

- div(A*(x,Lu)Vu*) = / in D, (3) 

with the same boundary conditions on u* as in (UJ). 

2.1.1. Random microstructure and deterministic homogenized 
coefficients 

Homogenization of elliptic equations with random coefficients has been 
extensively studied in the literature, and we refer to [251 [Ml 13 H3 for 
classical textbooks (see also the review article [5]). It is shown there that, 

if A e (x, u) — a(^ — ,u?J for some ergodic statistically homogeneous (i.e. 

stationary) random field A(x,uj) £ M. dxd (see e.g. [551111] for definitions), 
then the random solution u e (-,ui) to ([1]) converges, weakly in i/ 1 (D) and 
almost surely, to a deterministic function u* , solution to 

-div(A*Vu*) = / in D, 

with appropriate boundary conditions (say u* — on dD if (fT]) is comple- 
mented by u e (-,u) — on dD). The homogenized coefficient, denoted A* 
in the above equation, is a deterministic, constant matrix. 

In addition, the numerical procedure outlined above is a practical way 
to obtain a converging approximation of the homogenized matrix, in the 
sense that 

lim A*(x,uj) = A*, (4) 

almost surely, and for almost all x (see [TT]). Note that (@| can be equiva- 
lently written lim A*(x, uj) = A* for any fixed rj > 0. 

The only assumptions of ergodicity and stationarity do not allow for a 
precise convergence rate in ([4]) . If, in addition, one assumes that the matrix 



5 



A(x, uj) decorrelates at large distances at some given rate, then one can also 
obtain a convergence rate in ((4]) (see e.g. [28l [Tl] ) . A typical result is that 



E 



| A* (x, 



A* 



<C\~ 



a.e., 



(5) 



for some j3 > and C > that depend on the decorrelation rate, but are 
independent of x, rj and e, and where |-| is any norm on the dx d matrices. 

Note that, in the absence of ergodicity, the homogenized coefficients are 
a priori random matrices, that are invariant under the group of actions 
representing homogeneous statistical fields. 



2.1.2. Stochastic homogenized coefficients 

As we mentioned in the introduction, the Multi-Level Monte Carlo method 
is more efficient than a standard Monte Carlo method when the exact ho- 
mogenized coefficients are stochastic (otherwise, both methods are equally 
efficient). In stochastic homogenization, if no ergodicity is assumed, then 
the homogenized coefficients can be stochastic. In this work, we consider 
various cases in that setting. 

The first case we consider is when the coefficient in (TJ]) has the form 

A (x, — ,uj, uj'^J — A(x, uj) B ^— , uj'^j Id, 

where A and B are two random scalar valued functions and Id is the identity 
matrix. We thus see that uj corresponds to a randomness at the macroscopic 
scale, while uj' corresponds to a randomness at the microscopic scale. Let 
A* (x, lu, to') be the homogenized matrix, which depends on the macroscopic 
variables (x, uj), and also on the microscopic randomness uj' as no ergodicity 
is assumed on B. We will assume that 



A*(x,uj,uj')-H v (a (a 



X 

X, -,L),L)' 

e 



< C 



where the constant C and the rate /3 are independent of uj, x, e and rj. 

A second, more general case we consider is when the randomness does 
not explicitely split into a randomness at the macroscopic and the micro- 
scopic scales. The heretogeneous field in (TJ} then writes A (x, —,uj^j. We 

assume that A is scalar-valued, that we can do homogenization at every 
macroscopic point, and that the following assumption holds: 



E 



A*(x,uj) — H, 



(A (x, ■ 



< C 



for some constant C and rate /3 independent of x, e and ?y. This assump- 
tion is similar to the known results for ergodic homogeneous stochastic 
homogenization recalled in ([5]). 
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2.2. Multi-Level Monte Carlo approach 

We now briefly introduce the Multi-Level Monte Carlo (MLMC) approach 
in a general context. The reader familiar with this approach can directly 
proceed to Section HOI 

Let X(ui) be a random variable. We are interested in the efficient compu- 
tation of the expectation of X, denoted by E(X) . In our calculations below, 
X is a function of the homogenized coefficients or of the homogenized solu- 
tions. For example, we are interested in the expectation of the homogenized 
coefficients E(j4*), or in the two-point covariance function. In this case, 
we choose the random variable as X(u>) = [A*(xi,u)] i j [A* (x2,uj)] qp for 
some x\ £ D and X2 € D (and some components ij and qp of the homog- 
enized matrices). Other quantities of interest include e.g. statistics of the 
homogenized solution. 

To compute an approximation of E(X), a standard approach is the Monte 
Carlo (MC) method. One first calculates a number M of independent 
realizations of the random variable X (denoted X 1 , I < i < M) , and next 
approximates the expected value E(Jf) by the arithmetic mean (also called 
empirical estimator) : 



In this article, we are interested in Multi-Level Monte Carlo (MLMC) 
methods. The idea is to consider the quantity of interest X\ on different 
levels I. In our case, levels denote various representative volume sizes, or 
different mesh sizes. We assume that L is the level of interest, and that 
computing many realizations at this level is too computationally expensive. 
We introduce levels smaller than L, namely L — 1, . . . , 1, and assume that 
the lower the level is, the cheaper the computation of Xi is, and the less 
accurate Xi is with respect to Xl- Setting Xq = 0, we write 



The standard MC approach consists in working with M realizations of the 
random variable Xl at the level of interest L. In contrast, within the 
MLMC approach, we work with M; realizations of Xi at each level I, with 
Mi > M 2 > ■ ■ ■ > Ml- We write 



and next approximate E [Xi — Xi—i] by an empirical mean as above: 



Em(X) :=ttE^- 



L 



L 



E[X L ] =^E[I,-I,_ 1 ], 
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where X\ is zth realization of the random variable X computed at the level 
I (note that we have M[ copies of X\ and Xi-i, since M; < Mi-%). The 
MLMC approach consists in approximating E(Xl) by 



L 

E l (X l ):=J2Em, 
1=1 



(6) 



As will be seen below, the realizations of Xi used with those of Xi-\ to 
evaluate Em 1 (Xi — Xi-\) do not have to be independent of the realizations 
of Xi used with those of Xi + i to evaluate Em 1+1 — Xi) (see also 

Remark 15.11 below) . 

In the following, we are interested in the root mean square errors 



e-MLMc{Xi 



&MC {Xl 



E \\E(X L )-E L (X L ) 



E 



\&{X l )-E Ml {X l )\\ 



(7) 



(8) 



with an appropriate norm depending on the quantity of interest (e.g. the 
absolute value for any entry of the homogenized coefficient, the L 2 (D) norm 
for the homogenized solution). For the error estimation, we will use (see 
e.g. |13j ) that, for any random variable X, and any norm associated to a 
scalar product, 



E 



mx)~E M (x)\\ 



= — E 

M 



\X-E{X)\\ 



(9) 



2.3. Definition of meshes and representative volume sizes 

In our application, we will be dealing with various representative volume 
sizes, and also possibly various sizes of coarse meshes (see Figure [5] for 
illustration). In the framework of MLMC approaches, choosing a level 
I thus corresponds to choosing a particular RVE size, . . . We denote the 
hierarchy of coarse meshes on which we solve ([3]) by 

Hx>H 2 >-->H L . 

The number of realizations used at the level i for the coarse mesh size Hi 
is denoted Mj. We take 

Mi > M 2 > ■ ■ ■ > M L - 
As for the representative volumes, we take their sizes according to 

ill < m < ■ • ' < Vl 
and the corresponding number of realizations is denoted 

m\ > 77i2 ~> ■ ■ ■ > rriL- 
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One could also use various fine-scale meshes for solving the local represen- 
tative volume problems ([2]). We do not go in this direction in this work. 

Note that the level L always corresponds to the most expensive choice 
(large RVE, or fine mesh), and thus the smallest number of realizations. 
Note also that one does not have to take the same number of levels L for 
coarse-grid sizes and RVEs. 

H 

•p — coarse grid block 

^ffi RVE 

1 1 1 1 1 1 1 r ~~— — f me mesh 



Figure 2: Parameters in the numerical homogenization procedure. 



3. MLMC approach for the upscaled coefficients 

In this section, we describe how to use the MLMC approach to compute the 
upscaled coefficients defined in Section 12.11 and the two-point correlation 
functions. We focus on how to choose RVE sizes for the problems ©, and 
thus assume that these problems are exactly solved. Setting 



5 t (x) = JE \A*(x,-)-A*(x,-)\ 2 



where | • | is some matrix norm, we assume, following Section 12.1.21 that 

i /3/2 



Si(x)<C{-\ , (10) 



for some f3 > and C > independent of I, e, ?y and of the macroscopic 
point x £ D (in what follows, we keep the dependency with respect to x 
implicit in our notation). For some special cases, one can obtain an estimate 
for (3 rigorously. For more complicated cases, we suggest in Section 15.11 
below a pre-computation strategy that can provide an estimate for (3. Note 
that a Central Limit Theorem type result corresponds to (3 = d (see e.g. [5] 
for such estimates in a weakly stochastic case). 

For clarity, we summarize now our MLMC algorithm for the upscaled 
coefficients: 

1. Generate toi random variables u>\, . . . , tJ mi . 

2. For each level I, 1 < I < L, and each realization ujj, 1 < j < mi < mi, 
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Solve the RVE problems ([2]) on Y*: for any i = 1, . . . , d, 

dxv(A e (y,u>j)Vxi(y,Wj)) = in Y£, Xi{v,Uj) = Vi on dY*. 
Compute the homogenized matrix A*(x,u)j) with 

1 



VI < i < d, A\ (x,Uj)ei 



Vi Jv 



A e (y,uj J )Vxi(.y,u j )dy. 



3. For each level I, 1 < I < L, compute 



^ mi 



3=1 



4. Compute the MLMC approximation E L {A* L ) of the expected value 
E(A* L (x,-)) following ®: 



Let us now estimate the error in the approximation of E ^[-AJJ^.V for 

any entry ij (1 <i,j< d) of the matrix Aj?. To simplify the notation, we 
write the calculations below as if A£ were a scalar quantity independent 
of x. These calculations are to be understood as calculations on the entry 

For the MLMC approach, the error reads 



e M LMc(Al) = ,/E (E(AJ) - EHADY 



\ 



E 



.2=1 



1=1 



\ 



E 



L 

E 



<E\/ E (( E -- B rr.,)(A?-Af_ 1 ))' 
1 



E 



(a*-a*-i- e (a*-a-i))' 
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where we have used Q. Writing that A* - A*_ 1 = (A* -A*) + (A* -A*_ ± ), 
and since mi < m;_i, we deduce that 



&mlmc{A* l ) < 



1 



/mi 

1-1 

E 



E 



(A* - A* -E(A* - A*))" 



E 



^E^= 

L 9 

<E^= 



E 



(A* — A* ~E (A[ — A*)) 
1 



/mi 



(A* - A*) 
1 



/mi 



*l + -=a/E[(A*)2], 



/mi 



where, for ease of notation, we have introduced some mi+i < mj,. Us- 
ing ([ID]), we deduce that 



/3/2 



1 



/mi 



:VE P*) 2 



For a fixed error, the optimal choice for the number m,/ of realizations 
at level I (namely for the RVE of size r/i ) is reached when these error parts 
are equilibrated. Therefore, we choose 



mi 



) E[(A*) 




21 „ -2 , _ i 

1 7 ' — 



i > 



2 < / < £ + 1, 



(11) 



for some parameters a;, and we check that indeed tul+i < mi, provided 
ckl+1 = oljj. We then have 



/ £ x /3/2 i 
e M LMc(A* L ) < C — ) ^ a;. 



(12) 



z=i 



For comparison, we consider the error if we calculate the approximated 
upscaled coefficient only for the largest RVE (of size 771), using a standard 
MC method with mi independent samples. Using ©, we find that the 
MC error reads 



euc{A\ 



E 




(A* L — E (A* L )Y 



As pointed out above, A* is assumed to be a random quantity, with some 
positive variance. It is thus natural to assume that the variance of A* L is 
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roughly independent of L, and hence that the MC error is of the order of 
C/y/fhi,. To have an error of the same order as that given by the MLMC 

approach, we take rhj, = O ^ independent realizations. 

Now that we have chosen the number of realizations for both approaches 
so that they reach the same accuracy, we are in position to compare their 
cost. Let Ni denote the cost of solving the RVE problem ^ on the domain 
Y£ of size f]f. The number of degrees of freedom needed is of the order 
of {r)i/e) d . Assuming that N t = (mh) d , the MLMC cost is W^ MC = 



y mi Ni , hence 



1=1 

L 



w m L mc = E f^_\ af ^y + ^y E[iAr] a -> (2L) 



1=2 ^l-l 

In the case of the MC approach, the cost reads 

WMk MC 

On Figure [3l we plot the ratio — mc — f° r different numbers of levels L 

RVE 

and rates (3, with the choice rji = 2 l ~ L . Note then that the largest RVE 
is always of size t]l — 1, independently of L, and that the smallest RVE 
size depends on L, and is rji — 2 1_i . On the right plot, we consider the 
case when e is fixed at a very small value independent of L. This value 
is sufficiently small to ensure that, even for the largest considered L, the 
smallest RVE is larger than e (thereby ensuring scale separation). On the 
left plot, we consider a more practical situation (which is the regime we 
choose for our numerical experiments of Section [5J , when e depends on 
L and is always 10 times smaller that the smallest RVE. This leads to 
values of e that are larger (and thus easier to handle numerically) than 
that considered on the right plot. 

As we can see, for a given number L of levels, the larger the rate /3 is, the 
smaller the cost ratio is, at equal accuracy. Otherwise stated, the faster 
the convergence of the apparent homogenized matrix with respect to the 
RVE size, the more efficient the MLMC approach is. We also observe on 
the right plot that, at fixed /? and e, the gain in terms of cost first increases 
when L increases and then reaches a plateau for large L. 

Remark 3.1 

In the above calculations, we have assumed that the cost of solving a local 
problem scales linearly with the number Af of degrees of freedom. This is 
true if one uses iterative solvers and the condition number of the precondi- 
tioned system is independent of the small scale e. One can also compare the 
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cost between the MLMC and MC approaches under different assumptions 
(e.g. when the cost of solving a local problem scales as C(e)jV 1+7 for some 
7><U 

Remark 3.2 

To compute the optimal number of realizations following (|lip , one needs to 
know the value (3 of the rate in (110|) . In general, this rate is not analytically 
known. To address this difficulty, we propose in Section \5.1\ below some 
means to estimate the value of (3 based on a priori, offline computations. 

We have shown above how to estimate K(A* L (x, •)) at any macroscopic 
point x. Another important quantity is the two-point correlation function 

Cor A ,{x,y) :=E([A*(x,u)] tj [A*(y, U )]Jj 

between the components ij and qp of the homogenized matrix at points x 
and y (note that we work with non-centered values of A*). For simplicity, 
we only consider two fixed locations x and y. Consider mi independent 
realizations of the homogenized matrix A*^ at level I (1 < k < mi). We 
define 



1 r i r i 

Cor mi {A^):=—Y,W\x) A*> k {y) 

as an empirical estimator of E ( [Af(x, u)]. ■ [Af (y, u})} V The MLMC 
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approximation of the two-point correlation function Cor^* (x, y) then reads 



Cor L (A* L ) := £ (tfor m ,(AD - Cor^A^)) . 
i=i 

Remark 3.3 

We have considered above that we could exactly solve the RVE problems ([2]). 
In practice, these problems are solved numerically, within some accuracy. 
A natural extension of assumption (|10[) is to assume that the error in the 
approximation of A* ( due to working on a truncated domain Y ni of size r\i 
with a finite discretization on a mesh of size hj) satisfies 




for some constant C independent of hj, rji and e. An analysis similar to 
the one above then follows. Note also that it may be possible to solve the 
local problems on some RVEs with a coarser approximation and correct 
this using the nearby RVEs, computed at full accuracy, in the spirit of the 
strategy proposed in U2]j in another context. The adaptation of such an 
idea to our context goes beyond the scope of the current work. 

4. MLMC for the homogenized solution 

In this section, we show how to estimate the expectation of the homoge- 
nized solution using the MLMC approach. We also introduce an extension 
of that approach, namely the weighted MLMC approach, in Section |4~21 

4.1. Separable case 

In this section, we assume that the coefficient in (TTJ) reads 
Ae(x,u,u)') = A(x,u)B (-,&/) Id 

for two scalar valued functions A and B, and therefore satisfies a separation 
of scales assumption. The coarse-scale problem associated to the highly 
oscillatory problem (fTJ) is 



-div 



A(x,u) B* (uj')Vu* = f in D 



We expect most of the randomness of the coefficient at the coarse-scale to 
be in A(x, uS). We thus use a simplistic treatment for averaging over u/ and 
approximate the above equation by 

- div \l(x, ui) E u > [B*] Vu*j = / in D. (13) 



14 



We are going to compute an approximation of E(it*), using the tuples 
(Hi,Mi,rn,mi) for 1 < I < L. 
We first need to calculate the homogenized coefficient B*(lj'). To do so, 

we solve in each direction, 1 < j < d, and for each realization B % (— ,a/^ 
of the coefhcient, 1 < i < mi, the RVE problem 



-div 



B 



(V)v x ; 



in Y m , 



y on dY m , 



(14) 



and calculate the corresponding homogenized coefhcient: 

1 



[Bfl 



\Y, 



Note that we have kept implicit the dependency of Xj with respect to the 
level We then introduce 



mi f-f 

2 — 1 



which is an approximation (at level /) of 
introduce ui, solution to 



-div 



A(x,u)E mt (B?)Vui 



[B*]. We correspondingly 



/ in D. 



In turn, this equation is solved on a mesh of size Hi , for several realizations 
of A(x,u>). We thus eventually define uf (with 1 < / < L and 1 < k < Mi), 
solution (on a mesh of size Hi) to the coarse-scale equation 



div 



A\x,uj)E mi (BnVu l i 



(15) 



The expected value E(u/) is approximated in a standard Monte Carlo fash- 
ion by 

Mi 

E(ui) ^ E Ml ( Ul ) :=^$>*> 



fc=i 



where uf is the solution to (fT5l) . 

To approximate our quantity of interest, E(ul), we can hrst perform the 
above procedure only at the level L. This yields a standard Monte Carlo 
approximation of E(ui). 

An alternative approximation is that provided by the MLMC approach, 
which reads 



E l (ul) ■= y^-^M; (ui - ui-x) with m = 0. 
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Introducing the norm 



\X\\ = 



E[\\X(u)\\% HD) 



1/2 



the MLMC error is estimated following the same lines as in Section [3J We 
obtain 



\\E(u L ) - E L {u L )\\ <jr-^\\u*- Ul \\ +^ Wi \\u* 



1=1 



To bound from above \\u* — ui\\, we introduce u* Hi , approximate solution 
to (fT5|) on a mesh of size Hi. It follows that 

ii* ii^ll* * 1 1 i 1 1 * I 
\\u -ttj||<||u -u Hl \\ + \\u Hl -ui\\. 

The first term is a discretization error, which typically satisfies (e.g. if we 
use a PI Finite Element method) the bound \\u* — u* H || < Hi. For the 
second term, it holds (all expectations are taken w.r.t. u') 



i Hl - u i\ 



< 



< 



< 



E 



\E m (B{)-E(B*)\< 



E mi (B*)-E(B* 



\E(B?)-E(B*)\ 



\B* -E(B* 



Using our assumption PH|) . that is Sf 



for some f3 > 0, and 



assuming that the variance of Bf is essentially independent of I, we get 



\\E(u L )-E L (u L )\\<J2^( 



Hi+[- 



13/2 



c 



For the standard MC approach (with M independent samples), the error 
reads 

||E(umc) - E m( Um c)\\ < 



M 



provided the variance of umc is of order one. 



4.2. General, non-separable case 

In general, the coefficient in (JlJ is of the form A (x,ui, —J, where there is no 
separation between the macroscopic and the microscopic randomness. In 



1G 



this case, the RVE problems are parameterized by the macroscale position 
x, and thus need to be solved in each coarse-grid block (in contrast to the 
separable case considered in Section BTTl where the local RVE problem (fTH) 
is independent of x). 

At any level I, let N\ oc H^ d be the number of coarse-grid blocks. We 
denote by Vi the set of the macroscale grid points at which we solve a 
RVE problem, with Card Vi = JVj. We assume that the coarse grids are 
nested from one level to the other, so that V\ C Vi C • • ■ C Vl- As before, 
on each grid of size Hi, we solve Mi coarse grid problems. To calculate 
the effective coefficient, we solve the RVE problems at each coarse grid 
point and for each realization of A, and we next average the energy over 
the spatial domain. Since the sets Vi are nested, once we have computed 
A* Hi (using RVEs of size rji) at the macroscopic points of the coarse mesh 
of size Hi, we readily get A* H . for j < I (see Tabled]). Thus, at each 
level I < L, and at each point of the grid of mesh size Hi, we only have 
to solve Mi — Mi + i RVE problems (associated to independent realizations) 
on RVEs of size r\i, and not M; of them. 





H L 


Hl-i 


Hi 


# coefficients to calculate 
with RVE size rji 


m 






A* 


Mi - M 2 


VL-l 




A* 


A* 


M L -i ~ M L 


VL 


A* 


A* 


A* 


M l 


# coefficients 
on grid size Hi 


M l 


Ml- i 


Mi 





Table 1: Calculating the coefficients on the diagonal (shown in blue) will automatically 
give the lower triangular values in the matrix. 



We denote u*. H the solution to the coarse-scale equation discretized on 
a grid of size Hi , and where the effective coefficient is computed from local 
problems set on RVEs of size . 

To approximate E(u* £ Hl ), we can use a MLMC approach based on the 
solutions u*. H ., 1 < j < L. However, such an approach discards the 
solutions u*. H . for i < j, which are however easy to compute. Indeed, 
once the coefficient A*. H . has been obtained at some level j, computing 
the solutions u*. H . for all meshes i < j is as inexpensive as computing 
u r)j H on ly f° r the mesh j . 

To benefit from this fact, we can approximate E(u* H ) using a weighted 
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MLMC approach, which is defined as 

L 

Weighted ~ E "'^M, («? " («) 
i=l 

where a; (1 < I < L) are parameters to be determined, and 

j=i \ i=l 
1 L 

' 3=1 

where we have set M^+i = and u n .H = for any 1 < j < L. Note that 
if ai — 1 for all I, we recover the standard MLMC approach. 

Errors associated to the weighted MLMC approach are estimated in 
Appendix [21 



5. Numerical results 

We consider the problem 



-div 



.4 



^w^.u'jVue = / in D = (0, l) rf , 



complemented by boundary conditions that will be made precise below. 
Likewise, the function / will be given below. Note that the exact homog- 
enized coefficient is independent of these choices. 

In what follows, we compare our MLMC results with standard MC re- 
sults at the highest level. We equate the cost for calculating the coefficient 
and the solution separately and compare the errors (in contrast to the the- 
oretical analysis of Sections [3] and [H where we have equated the accuracies 
and compared the costs). 

We will consider both one-dimensional and two-dimensional examples. 
For the one-dimensional cases, we have implemented the method in Mat- 
lab, and used the analytical solutions of the various PDEs. In the two- 
dimensional cases, we use a rectangular mesh with cell-centered finite vol- 
umes. To solve the PDEs, we use the modular toolbox DUNE, the Dis- 
tributed and Unified Numerics Environment [51 M HH1 EH] • 

When we use the MLMC approach to approximate the homogenized 
coefficient, we consider L = 3 different RVE sizes r\i , unless specified oth- 
erwise. Likewise, when we compute the homogenized solution, we also use 
L = 3 different coarse grids of mesh size Hi. For all the computations, 
we have used the same fine grid (see Table [5]) . We have made sure that 
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I Hi hi r]i 



# cells in RVE of size r\i 




4096 



1024 



256 



Table 2: Parameters for the MLMC approach (two-dimensional cases). 

the smallest RVE we consider is much larger than the characteristic length 
scale e (given below for each example) of the field A. 

We explain in Section 15.11 how to numerically estimate the rate of con- 
vergence f3 in (fT0|) . In Section T5.21 we present numerical results for the 
homogenized coefficients. Next, in Section 15.31 we present numerical re- 
sults for homogenized solutions. 

5.1. Numerical study of the convergence rate 

In our theoretical study described above, we have assumed that 



for some constant C and rate f3 independent of e and r\ (see (|10l0 . In this 
section, we numerically estimate the parameter /3 on a practical example. 
The considered scalar coefficient A I — ,oj'j (defined for x £ D C K 2 ) is 



a random field with expected value E(A) = 10 (independent of x and e) 
and Gaussian covariance function 



with a = y/2, t = \/2 and r = er = 0.04 (recall that \x — x'\ denotes 
the Euclidean distance in R 2 ) . We generate samples of the coefficient with 
the Karhunen-Loeve expansion. By construction, the characteristic length 
scale e is related to the correlation length in cov(x, x'), which is of the order 
of er . 

For any 1 < I < L, we calculate the effective coefficients A*(u]j) for the 
RVE [0, r/i] 2 (with rji — 0.5 L ~') for various realizations 1 < j < mi- 



The theoretical reference value is A* — lim E (A*) , to which we cannot 
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access in practice. We thus define the reference value as 



1 J±. i "*L 



1=1 

where we have taken into account all the realizations on the RVEs [0, r/i} 2 , 
1 < I < L, in order to decrease as much as possible the statistical er- 
ror. In practice, we work with L = 4 and m = (mi, 11J2, m,3, TI4) = 
(2000,1000,300,140). For each RVE of size 1 < I < L, we expect 
from (Unj) that 



mi 



— XK« 

mi * — ' 



4* 



4* — 4* 



hence 



mi 



A 



mi 



re/ 1 



i=l 



/Sin - 



InC. 



(17) 



Results are shown on Figure |4j where we plot the computed data points 
(with error bars) and the corresponding linear regression line. We see that 
we find a straight line with slope (3 — 1.53 and intercept In C — 1.059 in the 
asymptotic regime 77 3> 1 . Note that the value of j3 is smaller than, but close 
to, the value /3 thco = d = 2 that would be obtained using a Central Limit 
theorem argument (see discussion below (flOl) '). In the numerical tests that 
follow, we will often consider only the three smallest RVE (rj = 0.5, 0.25 
and 0.125), for computational cost reasons. The slope of the regression line 
computed on the basis of these three smallest RVE decreases to j3 = 1.0095. 
These estimations will be useful in Sect ion f5 . 2 . 2 1 below (see Example 1). 

5.2. Computation of the homogenized coefficient 



We first consider the one-dimensional situation (Section I5.2.1[) and next 
turn to two-dimensional test cases in Section 15.2.21 

5.2.1. One dimensional examples 

Since the local problems (HJ are analytically solvable, we can afford to take 
many levels and many realizations at each level. 



Example 1 (separable coefficient) As a first test-case, we consider a 

ads 



coefficient A ( — uj' ) such that its inverse reads 



/V 



1 (f, W ,o/)= C + 5>(^)sin 2 - 
' ' L i=l ^ 



exp(w), 



20 




-4 -3.5 -3 , -2.5 -2 -1.5 



Figure 4: Using (|T7|) to estimate /?: computed data points (along with error bars) and 
the corresponding linear regression line with slope /3 = 1.53. 



where u> and Xi are i-i-d- random variables, uniformly distributed in [0, 1], 
ipi are fixed random numbers in [0.2, 2], and C > is a deterministic con- 
stant. Note that A -1 is uniformly bounded away from 0. This coefficient 
is separable in the sense that A -1 writes as a product of a function of u 
times a function of w' . For a fixed realization w, it is well known that, in 
the one-dimensional situation, the homogenized coefficient is the harmonic 
mean. Therefore the apparent homogenized coefficient on the RVE [a, b] is 




JV 



C(b-a) + Y,X*(u') 



lb — a sin (4nbtpi/e) — sin (ATratfi/e) 



{8iripi)/e 



0.5 L 

In our simulation, we use the values C = 1, N = 20 and e = (which 

ensures that the smallest RVE considered in the MLMC approach, of size 
0.5 L , is much larger than e, the characteristic length of the heterogeneous 
coefficient). As reference, we use the MC approach with 1000 realizations 
of the apparent coefficient on the largest RVE [ol,6l] = [0,0.5]. In what 
follows, a realization is determined by the tuple (oj, Xi(w')> • • • > Xjv (<**')) < 
Likewise, expectations are taken with respect to lo and lu' . 

For the MLMC approach, we use the RVEs [a h k] = [0, 0.5 L+1 "']. For 
this case, we expect that (3 = 2. Following (|TT|) . we hence take m = 
(4 l ~'to£, • • • ,4mt,fni) realizations. For comparison, we calculate the er- 
ror of the standard MC approach on the large RVE [<zl, br] = [0, 0.5], with 

772 b 

ffiL = — samples, so that both approaches share the same cost. 

bh 
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We are interested in comparing the relative mean square errors 



( e M l LMc) (At) 



(eMLMc(A* L )Y 



(e r M l c) (At) 



(eMc(Aj)) 2 

(mi)) 2 



where eMLMc(A* L ) and eMc(A* L ) are defined by ([7]) and ([8]). Since the 
errors depend on the set of chosen random numbers, we repeat the com- 
putations Nb = 10000 times and calculate the corresponding confidence 
intervals for the errors: 



mean[(e 



rel\2] 



1.96 std[(e 



rel \2l 



'Nb 



,mean[(e r ' e ')1 



1.96 std[(e re ') 2 ] 



'Nb 



We take L = 3, and show on Figure [5] the relative mean square errors on 
the expected value and the two-point correlation of the effective coefficient. 
For both quantities, we observe that the MLMC approach yields errors 2.5 
times smaller than the MC approach at equal computational work. 





15 20 25 30 35 



(a) Relative mean square errors on the expected (b) Relative mean square errors on the two-point 
value of the effective coefficient. correlation of the effective coefficient. 



Figure 5: Relative mean square errors with equated costs and m = ( 1 67713, 4m3, 7713) 
for the Example 1 (separable coefficient). 



Example 2 (separable stationary coefficient) We now consider an ex- 
ample where the effective coefficient docs not depend on a/, in the limit of 
infinitely large RVEs. We take A with inverse given by 

A^ 1 (x,uj,ui') = (C + ^2xi(u')l{i,i+i)(x) sin 2 (2nx) exp(w), 
V iez / 

where to and Xi are i-i.d. random variables, uniformly distributed in [0, 1], 
C = 1 and l[i ! i + i)(x) denotes the indicator function which is equal to 1 for 
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x E [i,i + 1) and to zero elsewhere. The apparent homogenized coefficient 
on the RVE [a, b] (to simplify, we choose a and b in Z) is 



A* j6 (w,w') = / A 1 {x,u,u J ')dx 



exp(w) 
b — a 



b-l 



C(6- a) +0.5 



In this case, the coefficient A is stationary in the variables (x,u)'), hence 
the standard stochastic homogenization theory holds: the exact effective 
coefficient is independent from uj' and reads 



A» = 



A 1 (x,ui,ui')dx 



[expH (C + 0.5E( X ))] 



Remark that, as expected, lim A* 6 (a;,a/) = almost surely in uj'. 

b — a— >oo ' 

In addition, the Central Limit Theorem holds for this case, thus /3 = 1 
in ([TOl). 

Following Sectional the theoretical reference value is E [A* L (uj, uj')} , where 
uj') is the apparent homogenized coefficient on the largest RVE. How- 
ever, this theoretical reference value is not easy to compute. We prefer to 
work with a different reference value, which is analytically computable, and 
which is very close to E [A* l (uj,uj')} when the RVE at level L is large. In 
the sequel, we use as reference 



A* 



E [A*(u)] = (C + 0.5E(x)) E [exp(-w)] 



1 - 1/e 
C + 0.25' 



By construction, A* e 



ref = lim E[Al(u,u')}. 

For the MLMC approach, we use the RVEs [a t ,bi] = [0, 100 x with 
m = (2 L ~'m£, • • • ,2mi,mi) realizations (recall that (3 — 1 in this case, 
and hence this choice for m agrees with (fTTj) '). Note that the smallest RVE 
is again much larger than the characteristic length scale of the field A. We 
compare this approach with a standard MC approach on the largest RVE 



[a L , b L ] = [0, 100 x 2 i " 1 ] that uses fh L 



approaches share the same cost). 

For this example, we have considered the choices L 



3, 5 or 7. On Fig- 



ure/ y 

-MC) 



ureO we compare the relative mean square errors {^mlmc) anc ^ ( e 
on the expected value and the two-point correlation of the effective coef- 
ficient (along with the corresponding confidence intervals obtained from 
Nb = 10000 different sets of random numbers). We again observe that the 
MLMC approach is more accurate (for the same amount of work), and that 
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the gain in accuracy increases if we increase the total number L of levels 
(this observation is consistent with Figure 13]). For L = 3, the gain is equal 
to 1.5 for both quantities, whereas it is equal to 3 for L = 5 and to 8 when 
L = 7. 

Remark 5.1 

On this example, we have also considered a MLMC approach where the 
realizations of Xi used in Emi{Xi — are independent from the real- 

izations of Xi used in Emi{Xi + \— Xi). More precisely (assuming L = 3 for 
the sake of simplicity) , this approach consists in approximating E(X^ = 3) 
by 

■ M 3 

2(^3(wi)-Xa(wi)) 
.i=i 

Mi 

M 1 -M 2 ^ V ' 

A ,:=i+m 2 

(18) 
rather than by 

M 3 

5^(X 3 (W<)-*2(W<)) 

.1=1 




as in We compare on Figure^ this method with a standard MC method, 
where the number of samples has been chosen to again equate the costs. We 
again observe that the MLMC approach (|18[) (with independent samples) 
is more accurate than the MC approach. We also observe that, at equal 
cost, a better accuracy is obtained when one uses (|19|) (with samples that 
are not necessarily independent) rather than (TTg)) . 



E ind{ x L=z) ■— 77- 



1 

M 2 - M 3 



M-2 



(X 2 (ui) - X^Ui)) 

U=l+M 3 



Example 3 (non separable coefficient) We now consider the coefficient 
defined by its inverse as 

A^ 1 (^-,u>,u)'^ = C(l +lo) + exp (uiu'sin ^-j) cos ^-J , 

where uj and w' are i.i.d. random variables uniformly distributed in [0.5, 1], 
0.5 L 

e = -jo - (the smallest RVE is thus large compared to e) and C = 2e 
(which ensures that A is uniformly bounded away from 0). This example 
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5 10 15 20 25 30 35 5 10 15 20 25 30 35 



(a) Relative mean square errors of the expected (b) Relative mean square errors of the two-point 
value of the effective coefficient (L = 3). correlation of the effective coefficient (L = 3). 




5 10 15 20 25 30 35 5 10 15 20 25 30 35 

m 5 m 5 



(c) Relative mean square errors of the expected (d) Relative mean square errors of the two-point 
value of the effective coefficient (L = 5). correlation of the effective coefficient (L = 5). 




5 10 15 20 25 30 35 5 10 15 20 25 30 35 

m 7 m 7 



(e) Relative mean square errors of the expected (f) Relative mean square errors of the two-point 
value of the effective coefficient (L = 7). correlation of the effective coefficient (L = 7). 



Figure 6: Relative mean square errors with equated costs and m = 
(2 L_l mi,'" ,2mi,mi), for the Example 2 (separable stationary 
coefficient). 
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5 10 15 20 25 30 35 



Figure 7: Relative mean square errors for the Example 2 (separable stationary co- 
efficient). The MLMC results have been computed following (Til?]) with 
(Mi, M2, M3) — (167713,47713,7713). Results following the approach described 
in Remark |5.1l labeled as 'MLMC ind', have been computed following (fl8|) 
with (Mi, M2, M3) = (I67713, 4777,3, m 7i)- The MC results have been computed 
with rfiMC samples so that the costs of 'MLMC, 'MLMC ind' and 'MC are 
equal. We work here with My_i = 4Mj (rather than Mj_i — 2Mj as in 
Figure [5]) to ensure that the number of samples per level decreases. 



is more challenging than the two previous ones as it is not separable. The 
apparent effective coefficient on [a, b] is 




As for the previous example, we use the practical reference value 

A* ref := BmE[^( U)U ')]=^. 

As for Example 1, we expect in this case that j3 = 2 and use the RVEs 
[<H, bi] — [0, 0.5 L+1 ~'] with m = (4 i_i 7nr J , ■ • • ,4t71l,777l) realizations for 
the MLMC approach, and compare its accuracy (at equal cost) with MC 
results on the RVE [ol, &l] = [0, 0.5]. Choosing L — 3, we show on Figure[5] 
the relative mean square errors (e r A f LMC ) and (e^| c ) on the expected 
value and the two-point correlation of the effective coefficient (confidence 
intervals have again been obtained from Nb = 10000 different sets of ran- 
dom numbers). Again, the MLMC approach yields an accuracy gain (here 
of the order of 2) over the MC approach, for both quantities. 
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Figure 8: Relative mean square errors with equated costs and m = (I6TO3, 4m3, TO3), 
for the Example 3 (non separable coefficient). 



5.2.2. Two dimensional examples 

We have seen in the previous section that the MLMC approach is efficient 
in the one dimensional case. We turn here to two dimensional test cases. 

Example 1 (separable coefficient) We first study the case when there is a 
separation in the randomness at the macroscopic level and the microscopic 
level. We set 

where A and B are both scalar valued. The random field B has expected 
value E(B) = 10 and a Gaussian covariance function: 

cov(x, x') = Cov 

with a = y/2, To = y/2 and r = ero = 0.04. We generate samples of the 
coefficient with the Karhunen-Loeve expansion. We take A(ui) = exp(aj), 
where uj is distributed according to the Gaussian law iV(0, 1). The ef- 
fective matrix is A*(u>,u>') = A(u>) B*(ui'). We only define levels I to 
approximate the expectation of B* . Thus, at each level I, we define 
Ai(lj,lj') = A(ui) B* (uj' '). Using mi independent samples at the micro- 
scopic level {^>j} 1< <m an d n x m i independent samples at the macro- 
scopic level ^ .„ , , we define, for any 1 < i < n, 

r L J J l<J<m; , l<i<n 1 ^ — — 3 

e^w) := — 5>?( w ><) = — S>K)5f(4). 

m .i — * j j Tii .i — * J J 
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To approximate expectations at the microscopic level, we use the MLMC 
approach, and introduce, for any 1 < i < n, 



L 



TO; 

i=l j=l 



Expectations at the macroscopic level are approximated using a standard 
MC approach on the macroscopic random variable oj. The reference quan- 
tity we are after is 



1=1 



which is in practice approximated by 

L 



11 1 ~ 

A ref = T E ~ E ^7 E A K) B *K0- 



Z=l 



In this case, the errors read 



'/ 



&mlmc{A* l ) — 



e M c{A* L 



\ 



1 ™ 



\ 



1 n 

-Y,[ A ref-En lL (Al)^) 
1=1 



where mi, is the number of samples used in the MC approach. As men- 
tioned above, we equate the computational work of the MLMC approach, 

mi y — J , with that of the MC approach, which is mi y — J . 



i=i 

This leads to taking 



Y.i=i m i faA) 2 

(We) 2 



We next compare the errors. We choose to work with L = 3 levels, 
n = 500 and, to compute the reference value A* e f, we used m ref = 

{m r 1 ef ,m r 2 ef ,m r 3 ef ) = (2000, 1000, 300). We also adopt the parameters 
of Table O On Figure |H we show the errors on the first entry of the effec- 
tive matrix, eAfc([^-I] n ) and e AfLMc([^I] n ), for m = (4to 3 ,2to 3 , m 3 ) 
(this choice is consistent with the value /3 = 1 in (fTU)k in turn, this as- 
sumption for j3 is consistent with our empirical estimation detailed in Sec- 
tion [5TTJ. We observe from these simulations that the MLMC approach 
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provides (roughly twice as) smaller errors than the standard MC approach 
for the same amount of computational work. Similar conclusions hold for 
the other entries of the effective matrix. 




Figure 9: Errors eMcd^Uii) and e MLMc{[A* L ] n ) for m = (4m 3 , 2m 3 , m 3 ), for 
Example 1 (separable coefficient). 



Example 2 (non-separable coefficient) We now consider a more difficult 
case, where there is no separation between uncertainties at the macro- and 
the microscopic levels. In general, such cases are difficult to handle, since 
having a sufficiently large number of samples to appropriately reduce the 
statistical noise is very expensive. We consider below a specific example for 
A such that we can solve the local problems (and thus compute the effective 
coefficient) analytically, due to the specific choice of boundary conditions 
in the local problem. Note that, in the limit of infinitely large RVEs, the 
effective coefficient does not depend on the precise choice of the boundary 
conditions set on the local problems (see [TTj'l. 
We consider the scalar coefficient (x = (xi,X2)) 

A (x,u, ^p^'J = M (xi,cj, -jp^'J M (x 2 ,uj, ' 

and write the local problems with Dirichlet and no-flow boundary condi- 
tions: 

-div (a(x,u),-, u'^j Vxi) 
n ■ Vxj 

with dY® = {x G dY n \xt — or x, = rj}. With these choices, the local 
problem reduces to a one-dimensional problem in the direction x 2 ; for the 



= OmY v = (p,rj) d , 
= Xi on 8Y n , 

= oonar^\ay D , 
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function \i that only depends on x\. For the first entry of the upscaled 
coefficient, we get 

A* n (uj,uj') = f~ J A^ 1 (xi,u, dx^j iy A 2 (x2,ui,^-,oj' 



dx 2 . 
(20) 



In our example we choose 

yl^f 1 ^w, — ,lu'^J = C(l +cj) +exp ^wo/sin ^— ^ cos ^— J , 
^4-2 ^2,W) = + exp(5w)) x 2 

+ (1 + x 2 )exp Hi + a; 2 )cjw'sin f— J J cos f— J , 

where u> and u/ are i.i.d. random variables uniformly distributed in [0.5, 1] 
and C = 2e. In this case, we see that (TITO holds with (3 — 2. To ensure 



scale separation even for the smallest RVE, we take e = — ■ 

To define the reference value of the effective coefficient, we run a MC 
approach on the RVE [0, 0.5] 2 with fh = 400000 realizations. It is possible 
to compute such a large number of samples in this two-dimensional test 
case thanks to the specific analytical expression ([2D) . 

The MLMC approach is run with L = 3 different levels, and m = 
(167713,47713,7713) realizations at each level (a choice which is consistent 
with (fTTj) and the fact that j3 = 2). To determine a confidence interval, we 
repeat the overall procedure with Nb = 2000 different sets of realizations. 
We compare on Figure [TU] the accuracies of the MC and MLMC approaches 
at equal computational cost. Again, the MLMC approach is more accurate, 
here by a factor roughly equal to 5. 

5.3. Computation of the homogenized solution 
5.3.1. One dimensional example 



As in Section 15.21 we start with the one dimensional situation where we 
know the reference solution exactly. To make the computations even sim- 
pler, we assume that, at the coarse-scale, the problem is subjected to ho- 
mogeneous Neumann boundary conditions. The coarse problem thus reads 

±(A*(x,u;,u>')^p\=f(x), K)'(O) = («*)'(!) = «*(Q) = 0, (21) 



dx \ dx r 

where the right-hand side satisfies / / = 0. The exact solution is 

Jo 

u*{x,uj,uj')= j (A*(t,u),u)')) F(t)dt, F(t)= f(z)dz. 
Jo Jo 
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Figure 10: Relative mean square errors with equated costs and m = (167713,47713,7713), 
for the Example 2 (non-separable coefficient). 



Let Xi denote the vertices of the grid, < i < N. The numerical ap- 
proximation of u* is a piecewise constant function, equal, on the interval 

(Xi-i,Xi), to 



= Y / (A*(x 3 ,u J ,u J ')r 1 f F{x)dx. 



In the spirit of the Example 3 in Section r5.2.1[ we assume that the apparent 
homogenized coefficient, obtained by solving the local RVE problem on 
[a, 6], reads 

{A*{x 1 lo,uj'))- 1 = C(l + exp(5w))x 



exp ( (1 + x)uico' sin I - I I — exp ( (1 + x)ojuj' sin I — 



e 

where w and uj 1 are i.i.d. random variables uniformly distributed in [0.5, 1] 
and C = 2e. We take f(x) = e x - (e - 1). 

The reference quantity is the expectation of the solution to (j2"Tj) . com- 
puted with the coefficient A*^ obtained by considering an infinitely large 
RVE: 

(^(x,^'))- 1 =C(l+eM^))x. 
This reference quantity reads 

E {K>o) = C f exp(x)a; - exp(a;) - (e - 1) + 1 

where C — C (l + 2 e^duA. 
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To compute an approximation of E(u^ ), we use the MLMC approach 
with L = 3 levels. The RVEs are defined by [a t ,bi] = [0, 0.5 L+1 ~ ; ] and the 
grid sizes are Sj — (0.25, 0.125, 0.0625). To ensure scale separation even for 

the smallest RVE, we take e = 

On Figure [TTJ the accuracy of the MLMC approach is compared with 
that of the MC approach at equal computational cost (error bars have been 
computed using Nb = 20000 different independent realizations of the whole 
computation) , for two choices of the number 9Jt of realizations at each level. 
We see that the choice 2Jt = (I6M3, 4M3, M3), which is consistent with the 
rate (3 = 2, yields the best results (and an accuracy gain of 33 %). 




M 3 M 3 



(a) 2rt = (I6M3, 4M 3 , M 3 ) (b) <M = (4M 3 , 2M 3 , M 3 ) 

Figure 11: Relative errors (in L 2 norm) of the solution (one dimensional example). 



5.3.2. Two dimensional example 

We now turn to an example in dimension two. The reference problem ([1]) 
is complemented with homogeneous (zero) Dirichlet boundary conditions, 
and the source term is f(x) = f(xi,X2) = 100(xi + £2). 

In the spirit of the Example 1 of Section 15.2.21 we take 
A(x,u,*,u/} = A(x,uj)B (-,<•</) . 

where A and B are scalar- valued, B is a log-normal distributed random 
field, B = e , with E(K) = and where the covariance function of K(x, ui') 

is cov(x,x') = a 2 exp (— l z ~f I j , with a = tq = y/2. The parameter e is 

such that €Tq = 0.04. The macroscopic random field is given by 

A(x,lu) =2 + \u>i sin(27rxi)| + |wa sin(27ra;2)| + \ijJ3 sin(7rxi)| 
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with independent and normally distributed Uk, 1 < fc < 3. 

Since the coefficient is separable, we are in the setting described in Sec- 
tion 14.11 In particular, the RVE problems are independent of the macro- 
scopic point x, and we can use the MLMC approach. For each level 
1 < I < L, we hence solve the coarse problem (fT5"]) on a grid of size Hi, for 
Mi realizations of A. This defines the solutions uf , 1 < k < Mi, 1 < I < L. 

The MC approach consists in working only at the level L, and thus 
solving, on a grid of size Hl, the problems 

M. 



-div [A k (x, w)£^(£|,)V<J = / in D, l<k< 

The reference solution is built as follows. At each level /, we first solve fj 15[) 
with mi — fh and 1 < k < Mi = M. The reference value is defined as 
the mean over both the levels and the number of realizations of all these 
solutions: 

L M 

i=i m k=i 

In practice, we take M — 1000 and fh = 50. 

We again work with L — 3 different levels and we equate the costs 
of the MC and the MLMC approaches for the computation of the ho- 
mogenized coefficients as well as that of the coarse scale solutions. This 
respectively implies that the parameters of the MC approach are fh = 
% 2 ( mi rjl + m 2 ri + m 3 7?|) and M = H% (M 3 H^ 2 + M 2 H^ 2 + M X H^ 2 ) . 

On Figure [T21 we show the relative L 2 -errors 



zmlmc{ul) 



eMc(u L ) 



\E^ l (ul)-E l (u l )\\ L 2 (d) 

)m f (u L )\\ LHD) 



\ E $ L (u L )-E^(u L )\\ LHD) 



Efl L {u L )\\ LHD) 



computed with the parameters 9K = {Mi,M2,M 3 ) = (32,32,16) and 
m = (rni,m 2 ,m3) = (50,40,20). Note that DJl is chosen based on the 
calculations presented in [1] (we have checked that these calculations also 
hold for finite volume methods). 

We actually repeat the whole procedure 200 times, and show on FigurelT2l 
the 200 values of the relative errors that we found. We see that these errors 
are essentially the same for all the realizations. A gain in accuracy of the 
order of 5 is obtained when using the MLMC approach, for an equal cost: 

E(e M LMc) w 0.1411, E(e M c) « 0.6851. 

The standard deviation of the MLMC error is also smaller: 

std M LMC = 0.0324, std M c = 0.0565. 



33 



0.9 



0.5 



0.8 



0.7 



0.6 




0.4 



0.3 




20 



40 



bC 



100 



120 



140 



I 60 



200 



Figure 12: Relative L 2 -errors Cmc and euhuc on the homogenized solution. We show 
the results for 200 different independent realizations. 
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A. Appendix: weighted MLMC approach 
analysis 



We estimate here the error associated to the weighted MLMC approach 
introduced in Section W% To this aim, it is useful to introduce the function 



L L 



i=i j=i 




E(m) 7^ K(u*). We successively estimate these two contributions. 
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Systematic error estimation Following the same lines as in Section [3l we 
obtain that 



\a-iu* - u\\ < y]£i,i [ai - aj+i] 



1=1 



where we have set q;l+i = and £j t i := \\u* — u*. H[ 



. Choosing now 



ai 



3=1 



to equilibrate the terms in the above error bound, we get 

L L 



\a\u — u\ 



< £l,l with «i = 2j 



£l,i 



i=i 



As shown in Section FTTl we have £jj < Hi + 



.7 = 1 



13/2 



■'3,3 



(22) 



For the standard MC approach, the systematic error reads 
\W-S^\\%H + S. 

To have the same systematic error, we choose the coarse grid size H 

L / \ P/2 L 



S y^ j {ai~ai + i)Hi andRVEs of size rj so that S = ( — J = y^(a; — ai + i)Si 



1=1 



i=i 



Statistical error estimation The statistical error of the weighted MLMC 
approximation satisfies 



1=1 



&1 

/Mi 



To equate the error terms in the above sum, we choose 

2 / 



Mi = C 



ai{Hi+S, 



for I > 2, Mi = C 



OJl 



7i(H + S) 



(23) 



for some constant C and some parameters 7;. The statistical error then 
satisfies 

11^ {E w * eighted ) — E w * eighted || = 0(H + S). 

For the MC approach, we choose M = C(H +S)~ 2 independent realizations 
and thus get a statistical error of the same order. 
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Cost comparison Now that we have chosen parameters such that the MC 
and the weighted MLMC approaches share the same accuracy, we are in 
position to compare their cost. 

As above, the cost of solving the coarse scale problems is 



wz~Z LMC = J2 M i H r 2 and w c ^; sc = mh- 



The dominating part of the computational cost however lies in solving the 
local RVE problems. For the MC approach, we assume that we need to 
solve these problems at N macroscopic points x. We thus have 

w£& = mn(*Y = c- ^ 



zj e 2 {H + 5f 



For the weighted MLMC approach, we assume that, at each level Z, we solve 
local RVE problems at N% < Hf 2 macroscopic points (with N% < N2 < 
■ ■ ■ < Nl). At each of these points, we only need to consider Mi — Mi + \ 
realizations. The computational work thus reads 



L 2 



N, 



1=1 



On Figure [TBI we show the ratio of the works for solving the coarse problems 
and the RVE problems, as a function of the number of levels L. The figure 
is made with the parameter ft = 2 (which corresponds to a Central Limit 
Theorem type convergence, see discussion below (JTUJ)) . As on Figure [31 
we consider two possible regimes for e. We see that a significant gain is 
achieved even for moderate values of L. 
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